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Abstract 

Within a simple model context, the sensitivity and stability of the ther- 
mohaline circulation to finite amplitude perturbations is studied. A new 
approach is used to tackle this nonlinear problem. The method is based on 
the computation of the so-called Conditional Nonlinear Optimal Pertur- 
bation (CNOP) which is a nonlinear generalization of the linear singular 
vector approach (LSV). It is shown that linearly stable thermohaline cir- 
culation states can become nonlinearly unstable and the properties of the 
perturbations with optimal nonlinear growth are determined. An asym- 
metric nonlinear response to perturbations exists with respect to the sign 
of finite amplitude freshwater perturbations, on both thermally dominated 
and salinity dominated thermohaline flows. This asymmetry is due to the 
nonlinear interaction of the perturbations through advective processes. 
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1 Introduction 



A recurrent theme in fundamental research on chmate variabihty is the sensitiv- 
ity of the ocean's thermohahne circulation. When state-of-the-art climate models 
are used to calculate projections of future climate states as a response to differ- 
ent emission scenarios of greenhouse g substantial spread in the model 

this spread is the diverse behavior of the 



results is found. One of t he reasons o : 

I 

thermohaline circulation (iMcAvaneyl . 



200l|). 



Stommel 



196ll ) where it is shown 



The sensitivity of the thermohaline circulation is caused by several feedbacks 
induced by the physical processes that determine the evolution of the thermoha- 
line flow. One of these feedbacks is the salt-advection feedback which is caused by 
the fact that salt is transported by the thermohaline flow, but in turn influences 
the density difference which drives this flow. The salt-advection feedback can be 
conceptually understood in a two-box model 
to cause multiple equilibria and hysteresis behavior. 

In many models of the global ocean circulation, it appears that several equilib- 
rium states may exist under similar forcing conditions. When the present equilib- 
rium state, with about 16 Sv Atlantic overturning, is subjected to a quasi-steady 
freshwater input in the North Atlantic, eventually the circulation may collapse. In 
this collapsed state, there is deepwater formation in the Southern Ocean instead of 
in the Nort h Atlantic and the fo r mation of North Atlantic Deep Water (NADW ) 



has ceased (jStocker et al. 



1992 



Rahmstorl . 



1995 



Manabe and Stouffer 



19991). 



As this multiple equilibrium regime seems to be present in many ocean models, 
it is important to determine whether transitions between the different states can 
occur due to finite amplitude perturbations. 

I n a variant of the Stommel-model for which the temperature relaxation is 



fast, 



Cessil (119941 ) studied the transition behavior between the different equilibria. 
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In this case, there are only three equihbrium states, of which one unstable. In the 
deterministic case, she finds that a finite amplitude perturbation of the freshwater 
fiux can shift the system into an alternate state and that the minimum amplitude 
depends on the duration of the disturbance. Regardless of the duration, however, 
the amplitude of the disturbance has to exceed a certain value for a transition 
to occur. Under stochastic white-noise forcing, there are occasional transitions 
from one equilibrium to another as well as fiuctuations around each state. 



In 



Timmermann and LohmannI (120001 ). the effect of multiplicative noise (through 



fast fiuctuations in the meridional th ermal tempe rature gradient) on the variabil- 



ity in a box model similar to that in 



Cessil (I1994J ) has been studied. It was found 



that the stability properties of the thermohaline circulation depend on the noise 
level. Red noise can introduce new equilibria that do not have a deterministic 
counterpart. 

Another line of studies uses box models that show intrinsic variability because 
of the existence of an oscillat o ry mo de in the eigenspectrum of the linear oper- 



ator. Griffies and TzipermanI (1l995l ) show that noise is able to excite an other- 



wise d amped eigenmode of the thermohaline circulation. 



Tziperman and loannou 



( 120021 ) study the non-normal growth of perturbations on the thermally driven 



state and identify two physical mechanisms associated with the transient ampli- 
fication of these perturbations. 

Stochastic noise can have a signi ficant effect on the raean states of the thermo- 



haline circulation and their stability 



2001 



Tziperman and loannoul . 



IHasselmannl . 



1976 



Palmer 



1995 



Velez-Belchi et al 



2OO2I ) . Some of these mechai iisms are intrinsically 



linear, such as the effects of non-normal growth considered in 



Tziperman and loannou 



( 120021 ). Others are essentially nonlin ear mechanisms, such as tho s e cau sing the 



noise-induced transitions reported in 



Timmermann and LohmannI (120001 ). 



To study linear amplification mechanisms, the linear singular vector (LSV) 
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metho d i s often used , with main apphcations to p r edicta bihty studies (IXue and Zebiak 



1997a 



b 



Thompson 



19981 ). iKnutti and Stockerl (120021 ). for example, found that 



the sensitivity of the ocean circulation to perturbations severely limits the pre- 
dictability of the future thermohaline circulation when approaching the bifurca- 
tion point. The LSV approach, however, cannot provide critical boundaries on 
finite amplitude stability of the thermohaline ocean circulation. 

In a system which potentially has multiple equilibria and internal oscillatory 
modes, its response to a finite amplitude perturbation on a particular steady 
state is a difficult nonlinear problem. In this paper, we determine the nonlinear 
stability boundaries of linearly stable thermohaline flow states within a simple 
box model of the thermohaline circulation. To compute these boundaries, we 
use the concept of the Conditional Nonlinear Optimal Perturbation (CNOP) and 
study optimal non/mear growth over a certain given time r. We extend results on 
linear optimal growth properties of perturbations on both thermally and salinity 
dominated thermohaline flows to the nonlinear case. We find that there is an 
asymmetric nonlinear response of these flows with respect to the sign of the 
finite amplitude freshwater perturbation and describe a physical mechanism that 
explains this asymmetry. 



2 Model and methodology 



a. Model 



To illustrate the a pproach, the the ory is applied to a 2-box model of the thermo- 



haline circulation (IStommel 



196ll ). This model consists of an equatorial box and 
a polar box which contain well mixed water of different temperatures and salin- 
ities due to an equatorial-to-pole gradient in atmospheric surface forcing. Flow 
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between the boxes is assumed proportional to the density difference between the 
boxes and, with a hnear equation of state, related to the temperature and salinity 
differences in the boxes. 

When the balances of heat and salt are nondimension alized, the gove rning 



Diikstral (|2000|)) can 



dimensionless equations (we use the notation in chapter 3 of 
be written as 

^ = Vi-ni+\T-S\) (la) 
^ = V2-Siv3+\T-S\) (lb) 

where T = Te — Tp, S = Se — Sp are the dimensionless temperature and salinity 
difference between the equatorial and polar box and = T — S* is the dimension- 
less flow rate. Three parameters appear in the equations ([T]): the parameter rji 
measures the strength of the thermal forcing, 772 that of the freshwater forcing and 
?73 is the ratio of the relaxation times of temperature and salinity to the surface 
forcing. Steady states of the equations ([T]) are indicated with a temperature of 
T, a salinity of S and a flow rate ^ = T — S . A steady state is called thermally- 
dominated when > 0, i.e. a negative equatorial-to-pole temperature gradient 
exists dominating the density. A steady state is called salinity-dominated when 
^' < 0, i.e. a negative equatorial-to-pole salinity gradient exists dominating the 
density. 

It is well known that the equations ([1]) have multiple steady states for certain 
parameter values. Here, we fix 771 = 3.0, 773 = 0.2 and use r/2 as control parameter. 
The bifurcation diagram for these parameter values is shown in Fig. [1] as a plot 
of ^ versus 772- Solid curves indicate linearly stable steady states, whereas the 
states on the dashed curve are unstable. There are thermally-driven (hereafter 
TH) stable steady states {'^ > 0) and salinity-driven (hereafter SA, ie., the 
circulation is salinity-dominated) stable steady states (^ < 0). The saddle- node 
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bifurcation points occur at 772 = 0.600 and ri2 = 1.052, and bound the interval 
in T]2 where multiple equilibria occur. Suppose that this bifurcation diagram 
represents both the present overturning state (on the stable branch with ^ > 0) 
and the collapsed state (on the stable branch with ^' < 0). To study the nonlinear 
transition behavior of the thermohaline flows from the TH state to the SA state 
and vice versa, we consider the evolution of finite amplitude perturbations on the 
stable states. 

The nonlinear equation governing the evolution of perturbations can be de- 
rived from equation ([l]). If the steady state {T,S) is given and T' = T — T, 
S' = S — S are the perturbations of temperature and salinity, then it is found 
that 



-(2|^| + 1)T' + sign{<I/)[T S' - ST' - T\T' - S') 



(2a) 



^ = -{2m+r]3)S' + stgn{^)[TS' - Sr - S'{T' - S')] (2b) 
at 

where sign{^) is sign of steady flow rate ^ . If the perturbations are sufficiently 
small, such that the nonlinear part of the equations ([2]) can be neglected, we find 
the tangent linear equation governing the evolution of small perturbations as 



dT - - - - 

'±L. = -(2m + l)T' + signimTS' -ST') 
at 

^ = -i2\^\+r^,)S' + stgnmfS'-ST') 
at 



(3a) 
(3b) 



b. Conditional Nonlinear Optimal Perturbation 



To study nonlinear mechanisms o. 
modified the LSV technique and 



amplification 



Mu 



Berloff and MeachamI (119961 ) 



(120001 ) proposed the concept of nonlin- 



ear singular vectors (NSVs) and nonlinear singular valu es ( NSVAs) 



cepts were successfully applied by 



Mu and Wangi (120011 ) and 



hese con- 



Durbianol (1200 ih to 



study finite amplitude stability of fiows in two-dimensional quasi-geostrophic and 



shallow- water models, respectively. In iMu and DuanI (120031 ). the concept of the 
conditional nonlinear optimal perturbation (CNOP) was introduced and applied 
to study the "spring predictability barrier" in the El Nino-Southern Oscillation 
(ENSO), using a simple equatorial ocean-atmosphere model. The "spring pre- 
dictability barrier" refers to the dramatical decline of the prediction skills for 
most of the ENSO models during the Northern Hemisphere (NH) springtime. 
The CNOP can also be employed to estimate the prediction errors of an El Nino 



or a La Nina event (IMu et al. 



20031). 



As readers may not be familiar with this concept, we give a brief introduction 
to CNOP. Considering the nonlinear evolution of initial perturbations governed 
by ([2]). In general, assume that the equations governing the evolution of pertur- 
bations can be written as: 

0, 



in 1] X [0, te] (4) 

a;|f=o = xq, 

where t is time, x{f) = (xi(t), X2(t), x„(t)) is the perturbation state vector 
and F is a nonlinear different iable operator. Furthermore, Xq is the initial per- 
turbation, X is the basic state, {x,t) G x [0,te] with Q a domain in i?", and 

te < +00. 

Suppose the initial value problem (jlj) is well-posed and the nonlinear propa- 
gator M is defined as the evolution operator of (jl]) which determines a trajectory 
from the initial time t = to time tg. Hence, for fixed te > 0, the solution 
x{te) = M{xo; x){te) is well-defined, i.e. 



X{te) = M{XQ-X){te) 



(5) 



So x{tf.) describes the evolution of the initial perturbation Xq. 

For a chosen norm || ■ || measuring a;, the perturbation Xq^ is called the 
Conditional Nonlinear Optimal Perturbation (CNOP) with constraint condition 
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I^Joll < , if and only if 



^(=^05) = max J(a;o) (6) 

1 1 X [) 1 1 < O 



where 

J{x,) = \\M{x,-x){Q\\ (7) 

The CNOP is the initial perturbation whose nonlinear evolution attains the 
maximal value of the functional J at time tg with the constraint conditions; 
in this sense we call it "optimal". The CNOP can be regarded as the most 
(nonlinearly) unstable initial perturbation superposed on the basic state. With 
the same constraint conditions, the larger the nonlinear evolution of the CNOP 
is, the more unstable the basic state is. In general, it is difficult to obtain an 
analytical expression of the CNOP. Instead we look for the numerical solution, 
by solving a constraint nonlinear optimization problem. 

To calculate the CNOP the norm 



\Xq\ 



is used. Using a fourth-order Runge-Kutta scheme with a time step dt = 0.001, 
perturbation solutions (T', S') are obtained numerically by integrating the model 
([2]) up to a time tg and the magnitude of the perturbation is calculated. Next, 
the Sequential Quadratic Programming (SQP) method is applied to obtain the 
CNOP numerically; the SQP method is briefly described in the Appendix. 

To compare the CNOPs with the LSVs, the latter are als o computed using 



the theory of linear singular vector analysis (I Chen et al.l . 119971 ). We also use the 
norm ([8]) in this analysis and first solve ([3]) to obtain the linear evolution of initial 
perturbations. Subsequently, the singular vector decomposition (SVD) is used to 
determine the linear singular vectors (LSVs) of the model. 
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3 Stability and sensitivity analysis 



In this section, we compute the CNOPs to study the sensitivity of the thermoha- 
hne circulation to finite amphtude freshwater perturbations in the two-box modeL 
Two problems are studied: (i) the nonlinear development of the finite amplitude 
perturbations for fixed parameters in the model and (ii) the nonlinear stability of 
the steady states as parameters are changed. Both thermally-dominated steady 
states {'^ > 0) and salinity-dominated ones {'^ < 0) are investigated. The ini- 
tial perturbation Xq is written as Xq = (Tq, S'q) = {6 cos 9 , 6 sin 9) , where 6 is 
magnitude of initial perturbation and 9 the angle of the initial vector with the 
X-axis. 

3.1 Finite amplitude evolution of the TH state 

For the thermally-driven stable steady state, we consider the state T = 1.875, S = 
1.275, ^ = 0.6 (shown as point "A" in Fig. [1]) with the fixed parameters rji = 3.0, 
?72 = 1.02, ?73 = 0.2. We choose = 2.5 and use 6 = 0.3 as a maximum norm 
(in the norm ([8])) of the perturbations. The time is about half the time the 
solution takes to equilibrate to steady state from a particular initial perturbation. 
The amplitude 6 = 0.3 is about 10% of the typical amplitude of the steady state 
of temperature and salinity (T, S). For 9 in the range 7r/4 < 9 < 5tt/A, the initial 
perturbation fiow has \1''(0) < 0. As this is typically caused by a freshwater fiux 
perturbation in the polar box, we refer to the perturbation as being of freshwater 
type. For other angles 9 G [0,27r] , the initial perturbation fiow has \E''(0) > 0, 
which is typically caused by a salt fiux perturbation in the polar box and we refer 
to it as being of salinity type. 

Using equation ([2]) and ([3]), both CNOPs and LSVs are computed versus 
the constraint condition 6, respectively. The numerical results, plotted in Fig. [21 
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indicate that the CNOPs are located at the circle ||a;o|| = S, which is the boundary 
of ball II^Coll < The directions of the LSVs, which are independent of 6, have 
constant values of 6i = 1.948 (dashed line) and 62 = 5.089 (not shown). The value 
of 6 for the CNOPs (solid curve) increases monotonically over the 6 interval 0.01 
to 0.3. The difference between CNOPs and LSVs is relatively small when 6 is 
small. 

Integrating the model ([2]) with CNOPs and LSVs as initial conditions, respec- 
tively, we obtain their evolutions at time tg, which are denoted as "CNOP-N" and 
"LSV-N"; these are shown in Fig. [3l For comparison, the linear evolution of the 
LSVs are also obtained by integrating the model ([3]) with the LSV as an initial 
condition; this is denoted as "LSV-L" in Fig. [31 It is clear that the evolution of 
the CNOPs is nonlinear in 6, while "LSV-L" only increases linearly. The line of 
"LSV-N" is between "CNOP-N" and "LSV-L", but the difference between "LSV- 
N" and "CNOP-N" is hardly distinguishable over the whole 6 interval. Though 
this difference is not significant in this TH state, it is significant in the following 
investigation for SA state. In fact, it is very hard to know this without previous 
calculation. 

Note that since our numerical results demonstrate that the CNOPs are all 
located on the boundary ||a;o|| = S, we are able to show the sensitivity of THC 
to finite amplitude perturbations of specific fixed amplitude 6 = 0.2. In Fig. H] 
the value of J at tg = 2.5 is shown for the linear and the nonlinear evolutions 
of the initial perturbations obtained by ([2]) and ([3]). For the linear case (dashed 
line in Fig. Hj), there are two optimal linear initial perturbations 61 = 1.948 and 
02 = 5.089 with the same value of J, J{5, 9i) = J (5, 62) = 0.16484, which are the 
LSVs. Note that 62 — 61 = it, which means that in the linear case perturbations 
with > and \I'' < behave similarly (and hence symmetrically with respect 
to the sign of \i/'). For the nonlinear model (solid curve in Fig. H]), there is one 
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global optimal nonlinear initial perturbation with 6*3 = 1.979, which is the CNOP, 
and one local optimal nonlinear initial perturbations at 9^ = 5.058, with values 
of 7(5,^3) = 0.22413 and J(5, ^4) = 0.13052, respectively. The results in Fig. 1 
for 6 = 0.2 coincide with the results in Fig. [2l 

There is another difference between the linear and nonlinear evolution of the 
perturbations. When the initial perturbations are freshwater (\E'' < 0), the non- 
linear evolution leads to a larger amplitude than the linear evolution. When the 
initial perturbations are saline (\E'' > 0), the nonlinear evolution leads to a smaller 
amplitude than the linear evolution. For example, the initial perturbations with 
di and ^3 are such that \E'' < 0, while the initial perturbations with 62 and 64 
have \E'' > 0. 

The values of J/ 6 obtained by integrating with the CNOPs as initial 
condition are shown for different S in Fig. [5|l. The corresponding evolution of 
is plotted in Fig. [Sb. To relate the result in Fig. 5a to previous ones, consider 
the value of J/(5 at t = 2.5 on the curve of 6 = 0.2. In Fig. 4, the maximum of 
J is 0.224 and hence J/6 = 0.224/0.2 = 1.12. It follows from the Figs. [5]that for 
the CNOP with 6 = 0.01, the flow rate recovers to the steady state ^ = 0.6 
rapidly. For the CNOP with a larger initial amplitude {6 = 0.1,0.2), it takes 
much longer for the thermohaline circulation to recover to steady state. This is 
different from a linea r analysis where the evo l ution is the same for all optimal 



initial perturbations (ITziperman and loannoul . 



2002). 



In summary, the results for the TH state show that fresh perturbations, with 
\E'' < 0, are more amplified though nonlinear mechanisms than saline perturba- 
tions, with \I'' > 0. This is consistent with the notion that perturbations which 
move the system towards a bifurcation point will be more amplified through 
non-linear mechanisms than perturbations that move the system away from a 
bifurcation point. 
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3.2 Finite amplitude stability of the SA state 

We consider the salinity-dominated SA state (f = 2.674, S = 2.796, ^ = -0.122) 
for a slightly smaller value of 772 than for the thermally-dominated TH state in 
the previous section (r/i = 3.0, ri2 = 0.9, 773 = 0.2). It is indicated as point "B" 
in Fig. [H Again for the time te = 2.5, using the corresponding equations ([2]) 
and (I3l), both CNOPs and LSVs are computed versus the constraint condition S, 
respectively and the results are plotted in Fig. [61 

Similar to the results in Fig. [2], the directions of the LSVs, which are indepen- 
dent of 6, have constant values of 61 = 2.796 and 62 = 5.938 (dashed line). The 
line for 61 is similar to 62, and is not drawn in Fig. [6l The direction of CNOPs 
(solid curve) increase monotonously with 6 varying from 0.01 to 0.125. Then, 
6 drops down to 2.857 at 5 = 0.125 and increase slightly with 6 in the interval 
0.125 < 6 < 0.17. After that, 6 jumps up to 5.195 aX 6 = 0.17 and increases 
slightly with S in 0.17 <S < 0.22. 

Using equations ([2]) and ([3]), the evolutions of both the CNOPs and LSVs 
of the SA state are shown in Fig. [71 All the three kinds of evolutions are ap- 
proximately the same when the initial perturbations are relatively small. But 
for larger 6 the difference between "LSV-N" and "CNOP-N" is remarkably larger 
than that of the TH state (Fig. [31). The values of J of "CNOP-N" are always 
larger than those of both "LSV-L" and "LSV-N" for 5 > 0.17. 

The values of J at a time tg = 2.5 for all 6 show (Fig. [8^) two optimal 
linear initial perturbations 61 = 2.796 and 62 = 5.938 with J{6,6i) = J(5, ^2) = 
0.0526. Again, because 62 — 61 = tt there is the symmetry in response with 
respect to the sign of \E''. For nonlinear evolutions, there is one globally optimal 
initial perturbation at 6^ = 5.246 (the CNOP) and two locally optimal initial 
perturbations at 64 = 0.251 and 9^ = 2.890, with J-values of J{6, 6*3) = 0.0963, 
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J(5, 6*4) = 0.0432 and J{5, 6^) = 0.0503, respectively. The initial perturbations 
with 61 and ^5 are of freshwater type (\E'' < 0), while the initial perturbations of 
02, O3 and 6'4 are of salinity type > 0). 

To understand the difference between the maxima located at ^3, 64 and 6^, a 
contour graph of J{6, 6) is drawn in Fig. [8)3. It is clear from Fig. [HJo that there 
are three groups of local maxima which are indicated by the dashed lines. When 
6 < 0.125 there are only two local maxima and the CNOP is located in the regime 
57r/4 < e < 2n. When 0.125 < 6 < 0.17, the CNOP jumps from the interval 
57r/4 < < 27r to the interval 37r/4 < ^ < 57r/4. This coincides with the jumping 
behavior of 6 in Fig. [6l When 6 > 0.17, there are three local maxima and the 
CNOP jumps from the interval 37r/4 < 6 < 57r/4 to a new interval 5 < 6 < 6. 
Both jumps are also shown in Fig. [6] and J has a remarkable increase after the 
second jump (Fig. [7]). 

Also for the SA state, the value of J/ 6 along trajectories obtained by integrat- 
ing equations ([2]) using the CNOPs as initial conditions are shown for different S 
in Fig. [9K. The corresponding evolution of is plotted in Fig. [9b. For t = 2.5 
and 6 = 0.2, the value of J/6 = 0.0963/0.2 = 0.48, where 0.0963 is the maximum 
in Fig. [8^. The flow rate recovers to the steady state (whose value is —0.122) 
shortly after being disturbed with a small amplitude CNOP {6 = 0.01). It takes 
much longer for the thermohaline circulation to recover to the steady state after 
being disturbed with a larger amplitude {6 = 0.1,0.2), respectively. The larger 
the CNOP is, the larger the transient effect. In contrast to the TH st ate, there 



i s now an oscillatory attraction to the SA state, already described in 



fll96lh 



Stommel 



In both the salinity-dominated SA state and the thermally-dominated TH 
state, the CNOP always moves the system towards the bifurcation point. The SA 
state (TH state) has an asymmetry in the nonlinear amplification of disturbances. 
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with a larger amplification for those with \E'' > (\E'' < 0). 

3.3 Sensitivity along the bifurcation diagram 

Even if a TH or SA state is linearly stable, it can become nonlinearly unstable 
due to finite amplitude perturbations. The methodology of CNOP provides a 
means to assess the nonlinear stability thresholds of the thermohaline flows; here 
this is shown for the two-box model (2). Thereto, we compute the CNOPs under 
different 6 constraints for linearly stable TH and SA states along the bifurcation 
diagram in Fig. [TJ 

Along the TH branch, we vary 772 from 1.043 to 1.046 with step 0.001 and 
thereby approach the saddle-node bifurcation (Fig. [1]). For each value of 772, the 
CNOPs are obtained under the constraint that the magnitude of initial pertur- 
bations is less 0.2 [6 = 0.2) and again the time tg = 2.5. The trajectories of 
the CNOPs are calculated by integrating the model ([2]) (Fig. [TOh). Next, the 
corresponding flow rates are drawn in Fig. [TOb . Both flgures indicate that the 
CNOPs damp after a while for the steady states labelled Ai, A2 and A^. While 
these three states are consequently nonlinearly stable, the CNOP for steady state 
^4 {V2 = 1.046) increases in time, which implies that this steady state is nonlin- 
early unstable (although it is linearly stable) to perturbations with 6 = 0.2. 

From the above results, it follows that for each value of 772, in the multiple 
equilibria regime of Fig. [H a critical value of 6, say 6c, must exists such that the 
TH state is nonlinearly unstable. 6c is deflned as the smallest magnitude of a 
flnite amplitude perturbation which induces a transition from the TH state to 
the SA state. The larger the value of 6c, the more stable the steady state is. 
Using the CNOP method, the values of 6c can be computed and the results for 
the TH states (from 772 = 0.95 up to the saddle- node bifurcation at 772 = 1.052 
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) are shown in Fig. [TTl The curve separates the plane into two parts. For the 
regime under the curve the steady state is nonhnearly stable and for the regime 
above the curve it is nonlinearly unstable. When the bifurcation point Q in Fig. [1] 
is approached, 5c decrease more and more quickly. The critical value 5c reduces 
to zero sharply, as 772 approaches the bifurcation point. This explains how the 
steady states lose their stability when the bifurcation point Q is reached. 

The same calculations are performed for steady states on the SA branch, when 
?72 varies from 0.75 to the value 0.60 at the saddle-node bifurcation. The CNOPs 
are obtained under the constraint that the magnitude of initial perturbations is 
less than 0.1 {5 = 0.1) with time te = 2.5. The trajectories of the CNOPs at 
each steady state (of which J and \E' are plotted in Fig. [T2|) indicate that the 
CNOPs damp for the steady states labelled Bi, B2 and B^. Although there is 
oscillatory behavior, these states are nonlinearly stable. However, the evolution 
of CNOP for steady state -B4 (772 = 0.70) increases, which implies that the SA 
state is nonlinearly unstable to this finite amplitude perturbation (although it is 
linearly stable). 

The critical boundaries 5c for the SA states (Fig. [T3l) show that 5c decreases 
monotonically with 772 and becomes zero at saddle-node bifurcation (772 = 0.6). 
Similar to Fig. [HI the curve separates the plane into two parts. For the regime 
under the curve the steady state is nonlinearly stable and for the regime above 
the curve it is nonlinearly unstable. When ri2 decreases from 1.0 to 0.6, the SA 
steady states approach the bifurcation point P in Fig. [H and the critical value 
5c tends to zero. This explains how the steady states lose their stability at the 
bifurcation point P. 
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4 Summary and Discussion 

Within a simple two-box model, we have addressed the sensitivity and nonlinear 
stability of (linearly stable) steady states of the thermohaline circulation. One of 
the remarkable results obtained by the CNOP approach is the asymmetry in the 
nonlinear amplification of perturbations with \E'' < (interpreted as a freshwater 
perturbation in the northern North Atlantic) and \1'' > (interpreted as a salt 
perturbation in the northern North Atlantic). 

When we use LSV analysis, there are two singular vectors xi and —xi that 
correspond to one singular value cxi (see Fig. [Hand Fig. [Hk)- If xi is a freshwater 
type perturbation ( \E'' < ), —Xi must be a salinity type perturbation ( \1'' > ). 
The conclusion from the linear analysis (using LSV) is that the thermohaline cir- 
culation is equally sensitive to either freshwater or salinity entering the northern 
North Atlantic. However, the nonlinear analysis (using CNOP) clearly reveals a 
difference in the response of the system to the two types of perturbations. 

The asymmetry can be understood by considering the nonlinear evolution of 
perturbations in the two-box model. For the TH state, according to ([2]), the flow 
rate perturbation v]/' = T' — S' satisfies 

^ = (25 - 2T - 1)T' + (2T -2S + r]s)S' - ^'^ (9) 

Integrating the above equation, we find 

^'(t) = ^'(0) + /* L(T', S')dT - f m'^dr (10) 

JQ Jo 

where L{T', S') is the linear part of ([9]). It is well known that the two line ar terms 
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in determine the linear stability of the steady state (iStommell . 

For an initial perturbation \E''(0) < (freshwater type), the nonlinear term is 
always negative and the freshwater perturbation is amplified. This is a positive 
feedback and the stronger the freshwater perturbation, the stronger the nonlinear 
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feedback destabilizing the TH steady state. A perturbation \l/'(0) > (salinity 
type) is damped by the negative definite nonhnear term. This is a negative 
feedback and the stronger the sahnity perturbation, the stronger the nonhnear 
feedback stabihzing the TH state. Nonhnear mechanisms hence make the TH 
steady state more stable to salinity perturbations. 

This explains the results in Fig. |H For the freshwater type initial perturba- 
tions (\E'' < 0), the nonlinear evolution of the initial perturbation of (j2]) is larger 
than the linear evolution of the initial perturbation of ([3]). For the salinity type 
initial perturbations (\1'' > 0), the nonlinear evolution is smaller than the linear 
evolution. In general, the nonlinear term yields positive (negative) feedback for 
negative (positive) \E'' in the case of thermally-dominated (TH) steady states. 

On the other hand, when the basic steady fiow is a SA state (^' < 0), then 
we have 

^ = {2f-2S- 1)T' + {2S -2f + 7]3)S' + (11) 
Similarly, we have, 

^'(t) = ^'(0) + /* L(r, S')dT + f ^'Ht (12) 
Jo Jq 

Hence, due to nonlinear effects the SA steady state becomes more unstable (sta- 
ble) to disturbances ^' > (\E'' < 0) which explains the results in Fig. ^1^. In 
general, the nonlinear term yields positive (negative) feedback for positive (neg- 
ative) \E'' in the case of salinity-dominated (SA) steady states. 

The physical mechanism behind the loss of sta bility of the TH state is often 
discussed in terms of the salt-advection feedback (jMarotzkd . 
ter (salt) perturbation in the northern North Atlantic decreases (increases) the 
northward circulation and hence decreases (increases) the northward salt trans- 
port. The salt-advection feedback is independent of the sign of the perturbation 
of the fiow rate \E''. In contrast, the nonlinear instability mechanism of the ther- 
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mohaline circulation depends on the sign of the perturbation of the flow rate 
as discussed above. 

The CNOP approach also allows us to determine the critical values of the finite 
amplitude perturbations (i.e. 5c) at which the nonlinearly induced transitions can 
occur. The techniques are currently being generalized to be able to apply them 
to models of the thermohaline circulation with more degrees of freedom. The aim 
is to tackle these problems eventually in global ocean circulation models. When 
applied to the latter models, the approach may provide quantitative bounds on 
perturbations of the present thermohaline flow such that nonlinear instability can 
occur. 
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Appendix: The SQP method 

The constrained nonlinear optimization problem considered in this paper, 
after discretization and proper transformation of the objective function F, can 
be written in the form 

min F(x), 

subject to 

Ci{x) < 0, for i = 1, 2, 3, • • • , n, 

where the Cj are constraint functions. It is assumed that first derivatives of the 
problem are known explicitly, i.e., at any point x, it is possible to compute the 
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gradient \/F(x) of F and the Jacobian J(x) = of the constraint 

functions q. 

The SQP method is an iterative method which solves a quadratic program- 
ming (QP) subproblem at each iteration and it involves outer and inner itera- 
tions. The outer iterations generate a sequence of iterates (x'^, A*^) that converges 
to (x*,A*), where x* and A* are respectively a constrained minimizer and the 
corresponding Lagrange multipliers. At each iterate, a QP subproblem is used to 
generate a search direction towards the next iterate A'^'^^). Solving such a 

subproblem is itself an iterative procedure, which is therefore regarded as a inner 
iterate of a SQP method. The following is an outline of the SQP algorithm used 
in this paper. 

Step 1. Given a starting iterate A°) and an initial approximate Hessian 
H^, set A; = 0. 

Step 2. Minor iteration. Solve by the following QP subproblem. 

min{[VF{x'')]^d'' + ^{d''^ H''d^)), 

subject to 

where d^ is a direction of descent for the objective function. 

Step 3. Check if [x^, A^) satisfies the convergence criterion, if not set x^'^'^ = 
x^ + ad^, where a <1. For A^"*"^, it is also deter mined by d^, which is automati- 



cally realized in the solver (IBarclay et al 



19971 ). Go to Step 4. 



Step 4. Update the Hessian Lagrangian by using the BFGS quasi-Newton 



formula (ILiu and Nocedal 



19891 ). Let = and y'^ = VL(x^+\ A'^+i) 



\/L{x'',X''), where VL = V-F -|- VcA. The new Hessian Lagrangian, H^^^, can 
be obtained by calculating 
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Then set k = k + 1 and go to Step 1. 

In the SQP algorithm, the definition of the QP Hessian L a grang ian is 
crucial to the success of an SQP solver. In iGill and Saunders! (1l997l ). H'' is a 
limited-memory quasi-Newton approximation to G = \/'^L, the Hessian of the 
modified Lagrangian. Another possibility is to define if^ as a positive-defi nite 



matrix related to a finite-difference approximation to G fjBarclay et al. 



19971 ). In 



this paper, we adopt the former one, which h as been shown to be e fficient for the 



nonlinearly constraint optimization problem (IMu and Duanl . 



20031 ). 
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Captions to the Figures 

Fig. 1 

The bifurcation diagram of the Stommel two-box model for rji = 3.0 and 773 — 0.2. 
The points labelled A and B represent the thermally-driven steady state and 
salinity-driven steady state, respectively, considered in section 3. The points la- 
belled P and Q represent the bifurcation points of the model, respectively. Solid 
curves indicate linearly stable steady states, whereas the states on the dashed 
curve are unstable. There are thermally-driven (TH) stable steady states > 0) 
and salinity-driven (SA, ie., the circulation is salinity-dominated) stable steady 
states (ff < 0). 
Fig. 2 

The values of 6 for both the linear singular vectors (LSV, dashed line) and for the 
Conditional Nonlinear Optimal Perturbation (CNOP, solid line) of the thermally- 
driven state under the conditions S < 0.3 and tg — 2.5. 
Fig. 3 

Values of J at the endpoints of trajectories at time = 2.5 for different values of 
S. These trajectories started either with Conditional Nonlinear Optimal Pertur- 
bation (CNOP-N, solid curve) and linear singular vectors (LSV-N, dash-dotted 
curve, hardly distinguishable from the solid curve) perturbing the thermally- 
driven state. Also included are the endpoints when the tangent linear model is 
integrated with the linear singular vectors as initial perturbation (LSV-L, dashed 
curve ). 
Fig. 4 

The magnitude of perturbations J obtained at te = 2.5 for the the evolutions of 
perturbations of the thermally-driven state in the tangent linear model and non- 
hnear model. The initial perturbations have the form (^'(0), S'{0)) — {S cos 6, 5 sin 9) 
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with 5 = 0.2. 
Fig. 5 

Values of (a) perturbation growth J/ 5 and (b) flow stream function ^ along the 
trajectories computed with Conditional Nonhnear Optimal Perturbation (CNOP, 
solid curve) initial conditions superposed on the thermally-driven state. The dif- 
ferent curves are for different values of S. 
Fig. 6 

The values of 9 for both the linear singular vectors (LSV, dashed line) and for 
the Conditional Nonlinear Optimal Perturbations (CNOP) of the salinity-driven 
state under the conditions S < 0.22 and te — 2.5 . 
Fig. 7 

The magnitude of perturbation J at the endpoints of trajectories at time tg — 2.5 
for different values of 6. These trajectories started either with Conditional Non- 
linear Optimal Perturbation (CNOP-N, solid curve) and linear singular vectors 
(LSV-N, dash-dotted line) perturbing the sahnity-driven state. Also included are 
the endpoints when the tangent linear model is integrated with the linear singular 
vectors as initial perturbation (LSV-L, dashed line). 
Fig. 8 

(a) Values of perturbation magnitude J obtained at tg = 2.5 for the the evo- 
lutions of perturbations of the salinity-driven state in the tangent linear model 
and nonlinear model. The initial perturbations have the form (T"(0), 5"(0)) = 
{5 cos 9 , 5 sin 9) with 5 = 0.2 . And (b) the contour plot of J{9,S), with contour 
interval of 0.02 from J = 0.02 to J = 0.08 (solid curve and dotted curve). The 
three group local maxima are indicated as the dashed curves . 
Fig. 9 

Values of (a) perturbation growth J/ S and (b) flow stream function ^ along the 
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trajectories computed with Conditional Nonlinear Optimal Perturbation (CNOP), 
initial conditions superposed on the salinity-driven state. The different curves are 
for different values of S. 
Fig. 10 

Values of (a) perturbation growth J and (b) flow stream function along the tra- 
jectories computed with Conditional Nonhnear Optimal Perturbation (CNOP), 
initial conditions superposed on the thermally-driven state for different values of 

Fig. 11 

The critical value of 6 {5c) versus the parameter controlling the thermally-driven 
state near the saddle-node bifurcation at r]2 = 1.05. 
Fig. 12 

Values of (a) J and (b) \E' along the trajectories computed with Conditional 
Nonlinear Optimal Perturbation (CNOP) initial conditions superposed on the 
salinity-driven state for different values of 772. 
Fig. 13 

The critical value of 5 {5c) versus the parameter controlling the salinity-driven 
state near the saddle-node bifurcation at 772 = 0.6. 
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Figure 1: The bifurcation diagram of the Stommel two-box model for rji = 3.0 and 
rjs = 0.2. The points labelled A and B represent the thermally-driven steady state and 
salinity-driven steady state, respectively, considered in section 3. The points labelled P 
and Q represent the bifurcation points of the model, respectively. Solid curves indicate 
linearly stable steady states, whereas the states on the dashed curve are unstable. There 
are thermally-driven (TH) stable steady states > 0) and salinity-driven (SA, ie., 
the circulation is salinity-dominated) stable steady states <Q). 
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Figure 2: The values of 9 for both the linear singular vectors (LSV, dashed line) and for 
the Conditional Nonlinear Optimal Perturbation (CNOP, solid line) of the thermally- 
driven state under the conditions S < 0.3 and te = 2.5. 
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Figure 3: Values of J at the endpoints of trajectories at time te = 2.5 for different 
values of 6. These trajectories started either with Conditional Nonlinear Optimal Per- 
turbation (CNOP-N, solid curve) and linear singular vectors (LSV-N, dash-dotted curve, 
hardly distinguishable from the solid curve) perturbing the thermally-driven state. Also 
included are the endpoints when the tangent linear model is integrated with the linear 
singular vectors as initial perturbation (LSV-L, dashed curve ). 
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Figure 4: The magnitude of perturbations J obtained at tg = 2.5 for the the evolutions 
of perturbations of the thermally- driven state in the tangent linear model and nonlinear 
model. The initial perturbations have the form (r'(0), S"(0)) = {6 cos 9,6 sm9) with 
6 = 0.2. 
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Figure 5: Values of (a) perturbation growth J/6 and (b) flow stream function along 
the trajectories computed with Conditional Nonlinear Optimal Perturbation (CNOP, 
solid curve) initial conditions superposed on the thermally-driven state. The different 
curves are for different values of 6. 
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Figure 6: The values of 9 for both the linear singular vectors (LSV, dashed line) and 
for the Conditional Nonlinear Optimal Perturbations ( CNOP) of the salinity-driven 



state under the conditions 5 < 0.22 and = 2.5 . 
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Figure 7: The magnitude of perturbation J at the endpoints of trajectories at time 
te = 2.5 for different values of 6. These trajectories started either with Conditional 
Nonlinear Optimal Perturbation (CNOP-N, solid curve) and linear singular vectors 
(LSV-N, dash-dotted line) perturbing the salinity-driven state. Also included are the 
endpoints when the tangent linear model is integrated with the linear singular vectors 
as initial perturbation (LSV-L, dashed line). 
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Figure 8: (a) Values of perturbation magnitude J obtained at tg = 2.5 for the the evo- 
lutions of perturbations of the salinity-driven state in the tangent linear model and non- 
linear model. The initial perturbations have the form (T'(0), S"(0)) = {dcos9,Ssm9) 

with 6 = 0.2 . And (b) the contour plot of J (9,6), with contour interval of 0.02 from 
J = 0.02 to J = 0.08 (solid curve and dotted curve). The three group local maxima are 
indicated as the dashed curves . 
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Figure 9: Values of (a) perturbation growth J/6 and (b) flow stream function along 
the trajectories computed with Conditional Nonlinear Optimal Perturbation (CNOP), 
initial conditions superposed on the salinity-driven state. The different curves are for 
different values of 5. 
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Figure 10: Values of (a) perturbation growth J and (h) flow stream function ^ along 
the trajectories computed with Conditional Nonlinear Optimal Perturbation (CNOP), 
initial conditions superposed on the thermally-driven state for different values of r]2- 
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Figure 11: The critical value of S (5c) versus the parameter controlling the thermally- 
driven state near the saddle-node bifurcation at 772 = 1.05. 
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Figure 12: Values of (a) J and (b) ^ along the trajectories computed with Conditional 
Nonlinear Optimal Perturbation ( CNOP) initial conditions superposed on the salinity- 
driven state for different values of 772 • 
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Figure 13: The critical value of S (5c) versus the parameter controlling the salinity- 
driven state near the saddle-node bifurcation at r]2 = 0.6. 
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